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Self-consistent simulations of star cluster formation from 
gas clouds under the influence of galaxy-scale tidal fields 
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ABSTRACT 

We present the first results of a project aimed at following the formation and long-term 
dynamical evolution of star clusters within the potential of a host galaxy. Here we focus 
on a model evolved within a simplified potential representing the Large Magellanic 
Cloud. This demonstrates for the first time the self-consistent formation of a bound 
star cluster from a giant molecular cloud. The model cluster reproduces the density 
profiles and structural characteristics of observed star clusters. 

Key words: stellar dynamics — methods: N-body simulations — globular clusters: 
general — open clusters and associations: general — galaxies: star clusters — stars: for- 
mation 



1 INTRODUCTION 

Star clusters are vital components of galaxies. Indeed, ex- 
tracting information from the massive globular clusters 
(GCs) of our Galaxy provides fundamental information on 
the epoch of galaxy formation. Furthermore, systems of ex- 
tragalactic GCs are used as key determinants for tracing 
the dynamical, chemical and gaseous evolution of their host 
galaxies (see Brodie & Strader 2006 for a review). The 
smaller open clusters also have a role to play as well, as 
understanding their destruction within in the Galactic disk 
impacts the growing field of Galactic archaeology (Freeman 
& Bland-Hawthorn 2002), for example. 

Efforts to model the evolution of individual star clusters 
have made great advances in the last decade. On the hard- 
ware front the advent of special-purpose GRAPE processors 
for calculating gravitational forces, with the latest incarna- 
tion being the GRAPE-6 (Makino 2002), have allowed N- 
body models of up to N ~ 100 000 stars to be completed in 
a reasonable timeframe (Baumgardt & Makino 2003). The 
next generation of GRAPE (available in 2008), the use of 
massively-parallel supercomputers, and even graphics pro- 
cessing technology, will push this N continually higher to- 
wards the realm of direct GC models. Complementing this 
is the push to make the models as realistic as possible by 
including algorithms to deal with processes such as stellar 
and binary evolution in concert with the treatment of the 
gravitational interactions of the stars. This is the case for 
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both Af-body (Hurley et al. 2005) and statistical (Fregeau 
& Rasio 2007; Giersz, Heggie & Hurley 2008) techniques. 

On a grander scale simulations of galaxy formation 
which identify the formation sites of star clusters have been 
performed (Bekki & Chiba 2005). Furthermore, the for- 
mation of stars within turbulent giant molecular clouds 
(GMCs) based at these formation sites can then be fol- 
lowed with smooth-particle hydrodynamics (SPH) simula- 
tions (Bekki & Chiba 2007). Our aim is to interface these 
galactic-scale simulations of star cluster formation with the 
latest modeling techniques for following the long-term evo- 
lution of star clusters. The immediate benefits will be: a) the 
first self-consistent simulations of star cluster evolution from 
formation through to destruction, and b) a non-simplistic 
picture of how GC systems evolve with time and how this 
impacts the interpretation of observed extragalactic GC sys- 
tems. This latter advance in particular will be important 
for understanding the connection between extragalactic GCs 
and galaxy formation by injecting much needed numeri- 
cal models in to what has become a data-dominated field 
(Brodie & Strader 2006). We will also be able to explore the 
origin of observed features of star cluster systems, such as 
the age and parameter distributions of the Large Magellanic 
Cloud (LMC) clusters (Mackey & Gilmore 2003). 

In this letter we provide an overview of the modeling 
process that will be used to achieve our goals. We illustrate 
this process by describing the evolution of a prototype model 
which demonstrates a working interface between the cluster 
formation and long-term evolution codes. It also entails the 
first simulation of bound cluster formation from a gas cloud 
within an external galactic potential. We then discuss future 
improvements to the model and details of how this fits in to 
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our plans for a full investigation of the long-term evolution 
of galactic and extragalactic star cluster systems. 



2 A TWO-FOLD MODELLING PROCEDURE 

The present model is two-fold as follows. We first deter- 
mine the spatial distributions of stars that are formed within 
GMCs gravitationally influenced by their host galaxies. In 
this first step, we investigate short-term star formation pro- 
cesses in GMCs based on GRAPE-SPH simulations of tur- 
bulent GMCs. Then we investigate the long-term dynamical 
evolution of collisional stellar systems, incorporating pro- 
cesses such as two-body relaxation and mass-loss from stars, 
for initial distributions of stars derived in the first step. In 
this second step, we emphasize that we are taking dynami- 
cally unrelaxed stellar systems formed from GMCs and plac- 
ing these within an iV-body code for following the evolution 
of very young SCs gravitationally influenced by the tidal 
field of their host galaxies. 

Previous work along these lines that bears mention 
starts with the simulations of galaxy formation reported by 
van Albada (1982, and see associated references within) and 
McGlynn (1984). These dealt with similar physical processes 

- the evolution of an initially irregular mass distribution 

- although the context and approach were quite different. 
Importantly these were collision-free simulations performed 
without the influence of an external potential. Aarseth, Lin 
& Papaloizou (1988) examined the collapse of proto-globular 
clusters through iV-body models. They took an initial distri- 
bution of fragments (representing low-mass pre-stellar sub- 
condensations) and documented the ensuing gravitational 
relaxation phase leading to the creation of a core-halo den- 
sity structure. These models also did not include an exter- 
nal potential. Furthermore, each of these studies did not 
model the fragmentation process leading to the formation 
of gas clumps or the subsequent creation of stars from these 
clumps. These prior studies do however provide a good basis 
to which we can compare the results of our models. 

2.1 Turbulent GMCs 

Radial mass density distributions, p(r), of GMCs with the 
sizes, r g , and masses, ra g , are assumed to have homogeneous 
spherical distributions. Here we set up the initial velocity 
fields due to turbulent flows within GMCs in the same way as 
Mac Low et al. (1998). We therefore assume that a turbulent 
velocity field within a CMC is a Gaussian random field with 
a power spectrum for ^ k < fc max as follows: 

P(k) = P k a , (1) 

where a is set to be 2.0 for most models and Po is a pa- 
rameter controlling the total kinematical energy due to the 
turbulent flow in the GMC. We mainly use fc max = 8 in our 
models, however, the final dynamical and chemical proper- 
ties of the simulated GCs do not depend strongly on fc max . 
The virial ratio, t v , is a free parameter described as: 

U = 2T k /W = /(Po) , (2) 

where Tk and W are the total kinematical energy and the 
absolute magnitude of the total potential energy for a GMC, 
respectively. As shown in equation (2), t v (0 < t v < 1) is 



determined by Po and thus can control the initial random 
kinematical energy of gas particles in the present study. We 
typically take t v = 0.25 in our models. Since an isothermal 
equation of state is suggested to be appropriate for star- 
forming interstellar clouds of molecular gas (e.g., Mac Low 
et al. 1998; Klessen, Heitsch & Mac Low 2000), we adopt 
the equation with the initial temperature of 10 K. 

2.2 Star formation and stellar feedback effects 

Gas particles with initial masses of rrii are assumed to be 
converted into new stellar particles if (i) local dynamical 
time scales are shorter than the local sound crossing time, 
and (ii) local gas densities exceed a threshold gas density, 
p th , of star formation (e.g., Nakasato, Mori & Nomoto 2000). 
Since the gravitational softening length, e, for gas particles 
of GMCs in a model is fixed during the evolution of the GMC 
(e ~ 10~ 2 r g ), there is a maximum density, p max , which the 
gas densities of individual SPH particles, pi, can not exceed 
in the adopted GRAPE-SPH method: p max ~ iV n m;/e g 3 , 
where N n is the total number of "neighbor particles" sur- 
rounding an i-th SPH particle. Here p max is estimated to be 
of the order of 10 3 atom cm -3 for models with m g ~ 1O 6 M . 
We thus assume that p t h = p max in the present study. It 
should be noted here that the above p max is significantly 
lower than the threshold gas density (~ 10 5 atom cm~ 3 ) of 
individual stars suggested by Elmegreen (2004). 

2.3 External gravitational potential 

As star clusters evolve they lose stars to their host galaxy at 
a rate that depends on the strength of the galactic potential 
and the orbit within this potential. Our plan is to improve 
the treatment of the external gravitational potential used 
in simulations of star cluster evolution by using the results 
of the galaxy-scale calculations to provide a 'live' model of 
the potential. However, to start with we adopt a smooth and 
static potential that is the current standard for iV-body star 
cluster models. We tailor this to an external gravitational 
potential reasonable for the LMC in our initial study as one 
of our goals is to compare our results with the observed 
properties of young star clusters in the LMC. 

Star-forming GMCs are thus assumed to be gravitation- 
ally influenced by the LMC represented by a point-mass, 
M ga i, of 0.9 x 1O 1O M . Within this simplified potential the 
GMCs are assumed to have circular velocities determined 
by their locations and the mass of the LMC. 

2.4 Dynamical evolution of SCs just after their 
birth 

The long-term dynamical evolution of the young SCs emerg- 
ing from the GMCs is then followed using the NB0DY4 code 
(Aarseth 1999, 2003). This code employs the fourth-order 
Hermite integration scheme, without softening, and is de- 
signed to exploit the fast evaluation of the gravitational 
force and its time derivative by the GRAPE-6. The use of 
the GRAPE-6 allows models of up to N = 100 000 to be 
completed although for N < 10 000 it is possible to make 
progress even when performing the full force calculation on 
the host computer. NB0DY4 includes algorithms for stellar 
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and binary evolution as described in Hurley et al. (2001). 
It allows for the full range of possible interactions within 
binary stars as well as dealing directly with the effects of 
close dynamical encounters: perturbations to binary orbits, 
collisions and mergers, formation of three- and four-body 
subsystems, exchange interactions, tidal capture and binary 
disruption. The tidal field of the host galaxy is modelled 
as described in Sec. 2.4 by placing the model cluster on a 
circular orbit at a specified radial distance, R sc , from the 
centre of a point-mass galaxy. As such NB0DY4 facilitates an 
investigation of the combined effects of internal (stellar, bi- 
nary and dynamical) and external (galactic tide) influences 
on the evolution history of a star cluster. 

The typical starting point for models of star cluster evo- 
lution is to assume that the star formation process is com- 
plete. An initial model is designed by first drawing the posi- 
tions of the stars from some density distribution: the n = 5 
polytrope proposed by Plummer (1911) and the family of 
King models (King 1962, 1966) derived from observations of 
mature star clusters are common choices, the former primar- 
ily because of its mathematical simplicity. Next the stellar 
velocities are determined by assuming the model cluster is 
in dynamical, or virial, equilibrium, i.e. 



h = V [W - 2E t ) = 0.5 



(3) 



where we have switched to the virial ratio definition used in 
star cluster simulations and Et represents the energy contri- 
bution from external forces (see Fukushige & Heggie 1995; 
Aarseth 2003, p. 11). The masses of the stars are assumed to 
be either equal or drawn from an initial mass function based 
on observations of field stars. In our new method we do away 
with these assumptions and feed the results of the star clus- 
ter formation calculations directly in to the NB0DY4 code. 
That is, the masses, positions and velocities of the newly 
formed stars. In this way we can ascertain if a bound cluster 
eventuates, while following the evolution of the density and 
velocity profiles of the young cluster. 



3 RESULTS OF A PROTOTYPE MODEL 

In this letter we highlight our procedure by describing the 
evolution of a proto-cluster of 8 420 equal-mass stars emerg- 
ing from a GMC. We set m g = 10 4 M Q and r g = 10 pc for 
the GMC and take the number of gas particles to be 2 x 10 4 . 
Thus the mass of the stars formed is 0.5M© per star. The 
timescale of star formation in this step is ~ 3 x 10 6 yr and the 
star formation efficiency is 42%. Full details of the formation 
process will be supplied in an upcoming paper. The spatial 
distribution of the proto-cluster is shown in the top panels 
of Figure Q] This is the initial model for the NB0DY4 step of 
the evolution and as such is given an age of OMyr. In this 
initial model there are no gas particles remaining but resid- 
ual gas can be communicated to the TV-body code, and its 
effect accounted for, in any future models. The stars in the 
initial model do not constitute a bound system (the total en- 
ergy is positive) and are definitely not in virial equilibrium: 
the ratio of kinetic to potential energy is 1.84 (this is Q v 
from eq. 3 with E t = 0). The half-mass radius, Th, is 5pc, 
the velocity dispersion is 2.8 km s" 1 , and the corresponding 
crossing-time for this initial model is approximately 7Myr. 
The cluster was placed on a circular orbit within our 
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Figure 1. Spatial distributions in the XZ- and XY-planes for the 
prototype model at ages of 0, 10, 30 and 300 Myr. The number of 
single stars, N B , and binaries, N^,, bound to the cluster at each 
age are denoted in the figure. 



simplified LMC potential at a distance of i? gc = 1.5 kpc from 
the galaxy centre. This matches the location of the parent 
GMC and corresponds to the scale length of the stellar disk 
of the LMC. The choice of tidal field gives an initial tidal 
radius of 8pc. We note that as the model cluster evolves, 
stars are denoted as having escaped from the cluster (and 
are removed from the simulation) when their distance from 
the cluster centre exceeds twice the tidal radius. This proto- 
type model was evolved to an age of 550 Myr using a single 
processor of the Swinburne supercomputer. At this point a 
bound cluster comprising 10% of the original stars remains. 

The total energy of the system first becomes negative 
after 5 Myr have elapsed. Figure [1] shows the spatial distri- 
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bution of the stars at ages of 10, 30 and 300 Myr. Already 
at 10 Myr the makings of a star cluster are evident and cer- 
tainly by 300 Myr the system has the regular appearance of 
a tightly bound star cluster. The virial ratio for the clus- 
ter stars (as defined by equation 3) is shown in Figure [2] 
for the first 100 Myr of evolution. We see that the ratio 
steadily decreases and that by an age of 30 Myr the clus- 
ter has reached approximate virial equilibrium. Certainly 
by an age of 55 Myr, when a two-body relaxation time (as 
measured at the half-mass radius) has passed, the cluster 
appears to be dynamically relaxed. When the model clus- 
ter is 100 Myr old the escape rate of stars from the cluster is 
14/Myr. This is down from an average of 80/Myr during the 
first 55 Myr of evolution but matches the escape rate from a 
comparison cluster started in virial equilibrium and evolved 
on the same LMC orbit. The model at 100 Myr contains 
2 817 bound stars. It has spatial parameters of r t = 5.6 pc, 
r\ — 1.2 pc and core-radius, r c = 0.36 pc. The half-mass re- 
laxation timescale has decreased to 28 Myr and the velocity 
dispersion is 1.5 km s . 

For comparison we have evolved the same initial model 
but at a distance of R sc = 8.5 kpc from the centre of a point- 
mass galaxy with M ga j = 9 x 1O 1O M0 - this resembles the 
orbit of an open cluster residing in the Solar neighborhood of 
the Galactic disk and is commonly referred to as a standard 
Galactic tidal field (Giersz & Heggie 1997) . Within a point- 
mass galaxy the tidal radius for a cluster of mass, M c , scales 
as 

r t oc (M c /M gal ) 1/3 fl gc . (4) 

So the Milky Way (MW) model with M ga i,MW = 
10Mgal,LMC and i? gc ,MW — 6i? gc ,LMC has r t ,MW — 3r t ,LMC 
and therefore represents a weaker tidal field. We see in Fig- 
ure [2] that this leads to the model approaching virial equi- 
librium on a similar timescale to that of the LMC model 
(20 — 30 Myr). To further understand the effect of the tidal 
field we have also evolved the initial model as an isolated 
cluster (Et = 0.0 at all times). This model approaches virial 
equilibrium on a timescale of ~ 40 Myr. Thus we see a weak 
trend for the virial timescale to increase for models with 
a decreased influence from the external tidal force. At the 
same time, the half-mass relaxation timescale is greater in 
models with a weaker tidal field, owing to larger N and rh 
at any particular time. We note that in equation (3) we have 
neglected a term that represents the rotation of the cluster 
(the Coriolis force) as this is negligible for clusters in dy- 
namical equilibrium (Fukushige & Heggie 1995). However, 
its absence explains the slight deviation from Q v = 0.5 for 
our relaxed models, as does the fact that a cluster can ap- 
proach a state of dynamical equilibrium but will never actu- 
ally reach it in practice. The Milky Way model at an age of 
100 Myr contains 6 724 stars and has r t = 21 pc, rh = 3.5 pc, 
r c = 0.36 pc and a velocity dispersion of 1.4kms~ . 

In Figure [3] we show the two-dimensional radial density 
profiles of the prototype cluster models (in the LMC tidal 
field) at 0, 30 and 300 Myr. The profiles are compared to 
two families of observationally determined density profiles: 
empirical King (1962) models based on Milky Way globu- 
lar clusters and Elson, Fall & Freeman (EFF: 1987) models 
based on young LMC clusters. The former are described by 
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Figure 2. Evolution of the virial ratio for stars in the prototype 
model (open circles). Also shown is the same model but evolved 
at a distance of 8.5 kpc from the centre of a point-mass Milky 
Way galaxy (red x symbols), and evolved as an isolated cluster 
with no external force (cyan open squares). 

a (f) = (TQ I -nr —p: I , (5) 

Ul + (r/rc) a ] Va [l + (rt/r c ) a ] V V 
and the latter by 

/ 2\-"t/ 2 

<r(r) = a (l + ^J , (6) 
where 

r c = a(2 2 /^l) 1/2 . (7) 

The main difference between the model profiles is the in- 
troduction of 7 in the EFF profiles which allows greater 
flexibility to fit the slope of the distribution exterior to the 
core for clusters that do not show significant tidal trunca- 
tion. As noted by Mackey & Gilmore (2003) the two families 
overlap if we set 7 = 2 in equation (6) and assume rt — > 00 
in equation (5). In fact, the profile of our prototype cluster 
at 30 Myr is best fitted by such a scenario: an EFF model 
with r c = 0.28 pc, 7 = 2 and a central surface density of 
2 500 stars pc -2 . The tidal radius at this time is 7pc but us- 
ing this in equation (4) gives a profile that is too truncated 
to fit the data for r > 1 pc. It is interesting to note that the 
functional form of the EFF profile was originally suggested 
by McGlynn (1984) to represent the equilibrium profile aris- 
ing from an initially irregular distribution of objects but in 
the context of dissipationless collapse of proto-galaxies. The 
profile at 300 Myr is taken as representative of the density 
profile at late stages in the evolution of the bound cluster. At 
more advanced times the profile becomes progressively nois- 
ier as the number of bound stars decreases. This profile is 
well fitted by an EFF model with parameters of r c = 0.16 pc, 
(Jo = 5 500 stars pc -2 and 7 = 2.6, noting that this is the me- 
dian slope determined by both EFF and Mackey & Gilmore 
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Figure 3. Radial surface density profile of the prototype LMC 
model at ages of (red x symbols), 30 (cyan circles) and 300 (open 
stars) Myr (projected along the Y-axis). Shown also are the best 
fitting EFF profiles for the 30 Myr data (cyan solid line: r c = 
0.28 pc, (j o = 2 500 stars pc" 2 , 7 = 2.0) and the 300 Myr data 
(black solid line: r c = 0.16pc, co = 5 500 stars pc -2 , 7 = 2.6) as 
well as the best fitting King (1962) model for the 300 Myr data 
(dotted line: r c = 0.16pc, cro = 5 500 stars pc -2 , rt = 5.0pc). 

(2003) from samples of young and intermediate age LMC 
star clusters. However, the data at 300 Myr are even better 
fitted by a King model with the same r c and o"o as well as the 
addition of r t = 5 pc (see Figure[3} . This is to be expected as 
the model cluster orbiting close to the centre of the LMC is 
strongly truncated by this stage. In comparison, the density 
profiles of the cluster evolved in the standard MW tidal field 
are adequately fitted across the 30 — 100 Myr timeframe by 
King models with r c — 0.26 pc and ao = 3 500 stars pc -2 . 

We note that the prototype LMC cluster reaches a min- 
imum in r c at an age of approximately 400 Myr which is as- 
sociated with the end of the core-collapse phase of evolution. 
The r c /rh value at this point is 0.06 which is a typical core- 
collapse value for standard models of star cluster evolution 
(see Hurley 2007). There are two binaries in the cluster at 
this age, both residing in the core. The maximum number 
of binaries in the simulation is recorded shortly after the 
cluster forms (5 Myr) when 36 binaries are present. 



its progenitor giant molecular cloud, approaches virial equi- 
librium on the order of a few crossing times and before a 
half-mass relaxation time has passed - a result that is free 
of the influence of the external tidal field. We have also found 
that the density profile of the model cluster resembles the 
distributions derived from observations - King (1962) and 
Elson, Fall & Freeman (1987) - from the time that virial 
equilibrium is reached and onwards. Idealized models that 
start in virial equilibrium can comfortably adopt these den- 
sity distributions when setting the initial positions of the 
stars. Our findings are distinct from previous simulations of 
the collapse of proto-globular clusters and proto-galaxies - 
which also observed the emergence of a core-halo structure 
at equilibrium - primarily in that the initial process of frag- 
mentation is explicitly followed using a three-dimensional 
hydrodynamic scheme and the influence of a galaxy-scale 
tidal field is modelled throughout the process. 

This letter describes a model that demonstrates a work- 
ing interface between simulations of star cluster forma- 
tion and long-term star cluster dynamical evolution. In this 
model we have included several aspects, such as the use of 
equal-mass stars and a point-mass galaxy, that are com- 
monly utilised in models of star cluster evolution but are 
not the most realistic approach. Future models will expand 
our study into a full investigation including: 

• the effect of a spectrum of stellar masses on the forma- 
tion and evolution of bound clusters - associated processes 
such as stellar evolution and mass segregation will affect the 
binding energy and the evolution timescales; 

• a variety of cluster sizes, in terms of physical size and 
also a greater number of stars; 

• an exploration of how the SFE assumed in the forma- 
tion step affects the final outcome, and also a greater under- 
standing of the conditions required for substantial binary 
formation; and, 

• a variety of CMC locations within the host galaxy and 
an extension to model the live gravitational potential of the 
host galaxy in concert with the internal star cluster evolu- 
tion, with dwarf, spiral and elliptical hosts considered. 

This will allow a much more detailed study of the effects of 
internal and external dynamical processes on the formation 
and evolution of star clusters in a galactic context. 
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4 SUMMARY 

We have demonstrated the self-consistent formation of a 
bound star cluster that exhibits all the markings of a reg- 
ular relaxed star cluster: dynamical equilibrium, core col- 
lapse, and an approximately spherical spatial distribution 
of the stars. This initially unrelaxed model, emerging from 
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